#load network
#load network
setwd("~/Desktop/Scent-Net")
gs4_auth(email=TRUE)
ℹ The googlesheets4 package is using a cached token for 'jbohey@uci.edu'.
net.long <- read_sheet("1plhxbVh5IQLNVMazqO4OvfM3X6vbll_Y6BTMJ94bkI4", na=c("","na", "NA")) %>%
mutate(pollinator_group= if_else(pollinator=="papilio_gothica", "butterfly", pollinator_group)) %>% #fix papilio_gothica from bee -> butterfly
mutate(plant= recode(plant,senecio_tall_int = "senecio_integerrimus", agoceris_glauca = "agoseris_glauca", agoceris_aurantiaca="agoseris_aurantiaca", hydrophyllum_fendeleri="hydrophyllum_fendleri"))
Auto-refreshing stale OAuth token.
✔ Reading from "caradonna_rmbl_interaction_networks_data_EDI".
✔ Range 'caradonna_rmbl_interaction_networks_data_EDI'.
net <- net.long %>% select(plant, pollinator, interactions) %>% pivot_wider(names_from = plant, values_from = interactions, values_fn= sum, values_fill= 0) %>% #pivots data set longer so plants as columns and pollinators as rows, replaces NAs with zeros
column_to_rownames("pollinator") #pollinators are in first column, makes them into rownames; #deletes the first pollinators column
net <- as.matrix(t(net)) # turns it from dataframe to matrix bc packages like matrix; t transposes (turns dataset sideways - rows -> columns; columns -> rows)
polls <- colnames(net) #now pollinators are column names
plants<-rownames(net)
snet <- sortweb(net) #sort by row & column totals; sorted rows by rowsums and columns by column sums. Biggest numbers in top left of matrix
#load scent
#Go to comm_vols for code to make scent.net
#Loads in vol, filters with bouquet, takes mean of each voc for each species
#From comm_vols.Rmd
scent.net <- read_csv("data/Volatiles_by_species_mean.csv") %>% column_to_rownames("species") %>% #species are now rownames
as.matrix() #network analysis likes matrix
Rows: 46 Columns: 356
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): species
dbl (355): .alfa.-Copaene, .alpha.-Cubebene, .alpha.-Farnesene, .beta.-Bisab...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
scent.net <- scent.net[,colnames(scent.net)!="Cyclohexane, 1,3,5-trimethyl-2-octadecyl-"] #filtering this compound bc it has 0 interactions or emissions... 355 vs 354. breaks code later
#loads in chemical class and gets rid of ?
compound.class <-read_sheet("1BlDIOs4STe5ZTTPYE6qez1Nsey_ovEpe5m-SN8Q55gg", sheet= "chemical_class") %>% mutate(major_class= factor(major_class,levels=c("aliphatic", "monoterpene", "sesquiterpene", "nitrogenous", "sulfur", "benzenoid")) %>% fct_na_value_to_level("other")) %>%
filter(compound_long!="Cyclohexane, 1,3,5-trimethyl-2-octadecyl-") #filtering this compound bc it has 0 interactions or emissions... 355 vs 354
! Using an auto-discovered, cached token.
To suppress this message, modify your code or options to clearly consent to
the use of a cached token.
See gargle's "Non-interactive auth" vignette for more details:
<https://gargle.r-lib.org/articles/non-interactive-auth.html>
ℹ The googlesheets4 package is using a cached token for 'jbohey@uci.edu'.
✔ Reading from "355 chemical compound class".
✔ Range ''chemical_class''.
#codes chemical class by colors
voc.class.colors <- setNames(brewer.pal(7,"Set3"), levels(compound.class$major_class))
##match paul’s plants with the plants we have vocs for
scentplants<-rownames(scent.net) #%>% tolower() %>% str_replace(" ", "_")
paulsplants <- plants
#plants not in both datasets
setdiff(scentplants, paulsplants) #gives us bonus plants
[1] "aquilegia_caerulea" "frasera_speciosa" "geum_aleppicum"
setdiff(paulsplants, scentplants)
[1] "hydrophyllum_fendleri" "phacelia_heterophylla"
bothplants <- intersect(scentplants, paulsplants) #plants found in amanda and paul data
polls.count <-count(net.long, pollinator_group,pollinator_family, pollinator) #counting number of rows for each pollinator
polls.class <- tibble(pollinator=polls) %>% left_join(polls.count) %>% pull(pollinator_group)
Joining with `by = join_by(pollinator)`
polls.class.dataframe <- tibble(pollinator=polls) %>% left_join(polls.count)
Joining with `by = join_by(pollinator)`
polls.family <- tibble(pollinator=polls) %>% left_join(polls.count) %>% pull(pollinator_family)
Joining with `by = join_by(pollinator)`
polls.order.all <- polls.class
#net.grp.all has all of pauls 45 plant species
#net.grp only has plant species found in Amanda's and paul's datasets
net.grp.all <- aggregate(.~ polls.class, as.data.frame(t(net)), sum) #lump network by pollinator class; plants columns, pollinators rows
rownames(net.grp.all) <- net.grp.all[,1] #make plants rows, pollinators column names
net.grp.all[,1] <- NULL
net.grp.all <- as.matrix(t(net.grp.all)) #transposes matrix
net.grp<- net.grp.all[bothplants,]
net.cut <- net[bothplants,]
##Pollinator groups
#Tells you plants major pollinator
plants.majorpoll <- setNames(factor(colnames(net.grp)[apply(net.grp, 1, which.max)]), rownames(net.grp))
#sets colors for ALL 5 pollinator groups
pcols.all <- setNames(brewer.pal(5,"Set2"), colnames(net.grp))
#Four colors for each major pollinator group
pcols <- pcols.all[ levels(plants.majorpoll)]
#gives % bee pollination for each plant species
plants.bees <- net.grp[,"bee"]/rowSums(net.grp)
plot(sort(plants.bees))
scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ] #cuts scent to match what paul has pollinator data for
# Cap of plants.bees
plant.bee.cap <- capscale(sqrt(scent.net.cut)~plants.bees)
anova(plant.bee.cap)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999
Model: capscale(formula = sqrt(scent.net.cut) ~ plants.bees)
Df Variance F Pr(>F)
Model 1 505636 1.0525 0.355
Residual 41 19696248
#p=0.348
plot(plant.bee.cap, type="text")
#gives % fly pollination for each plant species
plants.fly <- net.grp[,"fly"]/rowSums(net.grp)
plot(sort(plants.fly))
sort(plants.fly)
agoseris_aurantiaca amelanchier_alnifolia
0.000000000 0.000000000
castilleja_sulphurea gentiana_parryi
0.000000000 0.000000000
ipomopsis_aggregata lathyrus_lanszwertii
0.000000000 0.000000000
mahonia_repens mertensia_ciliata
0.000000000 0.000000000
ranunculus_inamoenus rosa_woodsii
0.000000000 0.000000000
vicia_americana viola_praemorsa
0.000000000 0.000000000
delphinium_nuttallianum hydrophyllum_capitatum
0.001774623 0.010204082
campanula_rotundifolia mertensia_fusiformis
0.011428571 0.011627907
erythronium_grandiflorum heliomeris_multiflora
0.019230769 0.033676658
sedum_lanceolatum heterotheca_villosa
0.035398230 0.044919786
lomatium_dissectum taraxacum_officinale
0.076923077 0.102841678
senecio_serra senecio_integerrimus
0.125000000 0.174044876
arenaria_congesta boechera_stricta
0.178947368 0.181818182
erigeron_speciosus solidago_multiradiata
0.202226345 0.205357143
claytonia_lanceolata helianthella_quinquenervis
0.244186047 0.267904509
fragaria_virginiana linum_lewisii
0.333333333 0.333333333
potentilla_hippiana erigeron_flagellaris
0.397425583 0.432900433
potentilla_gracilis eriogonum_umbellatum
0.433962264 0.442477876
pseudocymopterus_montanus draba_aurea
0.490196078 0.564285714
eriogonum_subalpinum androsace_septentrionalis
0.575000000 0.630000000
galium_boreale achillea_millefolium
0.875000000 0.942622951
agoseris_glauca
1.000000000
#gives % butterfly pollination for each plant species
plants.butterfly <- net.grp[,"butterfly"]/rowSums(net.grp)
plot(sort(plants.butterfly))
sort(plants.butterfly)
agoseris_glauca amelanchier_alnifolia
0.000000000 0.000000000
androsace_septentrionalis campanula_rotundifolia
0.000000000 0.000000000
castilleja_sulphurea draba_aurea
0.000000000 0.000000000
erythronium_grandiflorum fragaria_virginiana
0.000000000 0.000000000
gentiana_parryi hydrophyllum_capitatum
0.000000000 0.000000000
linum_lewisii lomatium_dissectum
0.000000000 0.000000000
mahonia_repens mertensia_ciliata
0.000000000 0.000000000
mertensia_fusiformis pseudocymopterus_montanus
0.000000000 0.000000000
ranunculus_inamoenus rosa_woodsii
0.000000000 0.000000000
sedum_lanceolatum senecio_serra
0.000000000 0.000000000
vicia_americana viola_praemorsa
0.000000000 0.000000000
heliomeris_multiflora potentilla_hippiana
0.003934730 0.006436042
delphinium_nuttallianum taraxacum_officinale
0.006654836 0.008119080
achillea_millefolium helianthella_quinquenervis
0.012295082 0.013262599
solidago_multiradiata potentilla_gracilis
0.013392857 0.018867925
heterotheca_villosa claytonia_lanceolata
0.022887701 0.023255814
lathyrus_lanszwertii arenaria_congesta
0.025000000 0.057142857
senecio_integerrimus erigeron_flagellaris
0.065494239 0.089466089
eriogonum_subalpinum galium_boreale
0.125000000 0.125000000
eriogonum_umbellatum erigeron_speciosus
0.150442478 0.179962894
ipomopsis_aggregata boechera_stricta
0.217391304 0.363636364
agoseris_aurantiaca
1.000000000
#gives % hummingbird pollination for each plant species
plants.hummingbird <- net.grp[,"hummingbird"]/rowSums(net.grp)
plot(sort(plants.hummingbird))
sort(plants.hummingbird)
achillea_millefolium agoseris_aurantiaca
0.000000000 0.000000000
agoseris_glauca amelanchier_alnifolia
0.000000000 0.000000000
androsace_septentrionalis arenaria_congesta
0.000000000 0.000000000
campanula_rotundifolia castilleja_sulphurea
0.000000000 0.000000000
claytonia_lanceolata draba_aurea
0.000000000 0.000000000
erigeron_flagellaris erigeron_speciosus
0.000000000 0.000000000
eriogonum_subalpinum eriogonum_umbellatum
0.000000000 0.000000000
erythronium_grandiflorum fragaria_virginiana
0.000000000 0.000000000
galium_boreale gentiana_parryi
0.000000000 0.000000000
helianthella_quinquenervis heliomeris_multiflora
0.000000000 0.000000000
heterotheca_villosa hydrophyllum_capitatum
0.000000000 0.000000000
lathyrus_lanszwertii linum_lewisii
0.000000000 0.000000000
lomatium_dissectum mahonia_repens
0.000000000 0.000000000
mertensia_ciliata mertensia_fusiformis
0.000000000 0.000000000
potentilla_gracilis potentilla_hippiana
0.000000000 0.000000000
pseudocymopterus_montanus ranunculus_inamoenus
0.000000000 0.000000000
rosa_woodsii sedum_lanceolatum
0.000000000 0.000000000
senecio_serra solidago_multiradiata
0.000000000 0.000000000
taraxacum_officinale vicia_americana
0.000000000 0.000000000
viola_praemorsa senecio_integerrimus
0.000000000 0.002425713
boechera_stricta delphinium_nuttallianum
0.090909091 0.224933452
ipomopsis_aggregata
0.608695652
#gives % wasp pollination for each plant species
plants.wasp <- net.grp[,"wasp"]/rowSums(net.grp)
plot(sort(plants.wasp))
sort(plants.wasp)
achillea_millefolium agoseris_aurantiaca
0.0000000000 0.0000000000
agoseris_glauca amelanchier_alnifolia
0.0000000000 0.0000000000
androsace_septentrionalis arenaria_congesta
0.0000000000 0.0000000000
boechera_stricta campanula_rotundifolia
0.0000000000 0.0000000000
castilleja_sulphurea claytonia_lanceolata
0.0000000000 0.0000000000
delphinium_nuttallianum draba_aurea
0.0000000000 0.0000000000
erigeron_flagellaris erigeron_speciosus
0.0000000000 0.0000000000
eriogonum_umbellatum erythronium_grandiflorum
0.0000000000 0.0000000000
fragaria_virginiana galium_boreale
0.0000000000 0.0000000000
gentiana_parryi helianthella_quinquenervis
0.0000000000 0.0000000000
heliomeris_multiflora hydrophyllum_capitatum
0.0000000000 0.0000000000
ipomopsis_aggregata lathyrus_lanszwertii
0.0000000000 0.0000000000
linum_lewisii lomatium_dissectum
0.0000000000 0.0000000000
mahonia_repens mertensia_ciliata
0.0000000000 0.0000000000
mertensia_fusiformis potentilla_hippiana
0.0000000000 0.0000000000
pseudocymopterus_montanus ranunculus_inamoenus
0.0000000000 0.0000000000
rosa_woodsii sedum_lanceolatum
0.0000000000 0.0000000000
senecio_integerrimus senecio_serra
0.0000000000 0.0000000000
taraxacum_officinale vicia_americana
0.0000000000 0.0000000000
viola_praemorsa heterotheca_villosa
0.0000000000 0.0005347594
solidago_multiradiata potentilla_gracilis
0.0178571429 0.0377358491
eriogonum_subalpinum
0.0750000000
##Heatmap1
#pp_heatmap
heatmaply(net, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
row_side_palette=function(n) pcols,
col_side_colors=data.frame(PollinatorGroup=polls.class),
col_side_palette=function(n) pcols.all)
Warning in doTryCatch(return(expr), name, parentenv, handler): unable to load shared object '/Library/Frameworks/R.framework/Resources/modules//R_X11.so':
dlopen(/Library/Frameworks/R.framework/Resources/modules//R_X11.so, 0x0006): Library not loaded: /opt/X11/lib/libSM.6.dylib
Referenced from: <FFA47D77-8F35-36FC-B0E5-38351B8D9512> /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/modules/R_X11.so
Reason: tried: '/opt/X11/lib/libSM.6.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/opt/X11/lib/libSM.6.dylib' (no such file), '/opt/X11/lib/libSM.6.dylib' (no such file), '/Library/Frameworks/R.framework/Resources/lib/libSM.6.dylib' (no such file), '/Library/Java/JavaVirtualMachines/jdk-11.0.18+10/Contents/Home/lib/server/libSM.6.dylib' (no such file)
##Matrix plot
#interaction matrix plot
visweb(net, type="nested", circles=T, boxes=F, circle.max=1.2)
##Web plot
#fix names of species (plants rows, polls cols)
net.nice.labels <- net
rownames(net.nice.labels) <- str_to_sentence(str_replace_all(rownames(net), "_"," "))
colnames(net.nice.labels) <- str_to_sentence(str_replace_all(colnames(net), "_"," "))
#interaction web plot
par(mar=c(0,0,0,0))
plant.pollinator.network <- plotweb(sqrt(net.nice.labels), empty=F, text.rot=90, method="cca", col.high=pcols.all[polls.order.all], col.low=pcols[plants.majorpoll], col.interaction=pcols.all[polls.order.all], bor.col.interaction=pcols.all[polls.order.all], bor.col.high=pcols.all[polls.order.all], bor.col.low=pcols[plants.majorpoll])
library(ggplot2)
ggsave(plant.pollinator.network, file='plant.pollinator.networkplot.jpg', width=130, height=105, units= 'mm', dpi=300)
#Network topography
#computes network statistics (like nestedness)
#networkstats <- networklevel(net)
#write_csv(enframe(networkstats), "data/p-p network stats.csv")
#computes network statistics for voc-pollinator network
#scentnetworkstats<- networklevel(scentpollnet)
#write_csv(enframe(scentnetworkstats), "data/scent network stats.csv")
scentnet.stats <- read_csv("data/scent network stats.csv")
Rows: 48 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): name
dbl (1): value
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#the network statistics null model: TODO- broken
#network.nullmodel<-nullmodel(scentpollnet)
#nullmodel.stats<- map(network.nullmodel, ~networklevel(.x )) #network.nullmodel[1:2]; index=c("NODF", "weighted NODF", "connectance", "H2", "modularity Q")
#nullmodel.stats %>% map_dfr(~enframe(.x), .id = "run") %>% write_csv("data/scent network null model statistics.csv")
nullmodel.stats <- read_csv("data/scent network null model statistics.csv")
Rows: 48000 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): name
dbl (2): run, value
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#idea - pivot wider so network properties are the columns, rows are run number, then average by values in column -> pivot longer or store in new variable
#nullmodel.stats %>% pivot_wider(names_from = "Name", values_from = "Value") %>% column_to_rownames("stats") %>% as.matrix()
null.NODF <- read_csv("data/scent network null model statistics.csv") %>% filter(name=="NODF")
Rows: 48000 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): name
dbl (2): run, value
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ggplot(null.NODF, aes(x=value))+ geom_histogram() + geom_vline(aes(xintercept=filter(scentnet.stats, name=="NODF") %>% pull(value)))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
##Heatmap2
scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ] #cuts scent to match what paul has pollinator data for
plants.majorpoll.cut <- plants.majorpoll[rownames(scent.net.cut)] #cuts Pauls data to match Amanda's
plants.bees.cut<- plants.bees[rownames(scent.net.cut)]
heatmaply(sqrt(scent.net.cut), scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
row_side_colors = data.frame(major.polls=plants.majorpoll.cut),
row_side_palette=function(n) pcols )
Warning in distfun(x): you have empty rows: their dissimilarities may be
meaningless in method "bray"
#col_side_colors=data.frame(PollinatorGroup=scent.net), #put compound class here
# col_side_palette=function(n) pcols.all)
#Distance Matrix (for bray curtis distance) #Mantel test
#scent-plant network
distance.matrix <- vegdist(scent.net.cut)
View(as.matrix(distance.matrix)) #tells you how different the smells are. If 0 there is no difference, scents are the same; 1 = no compounds in common
#plant-pollinator network
#how similar are the plants in terms of their pollinator array/composition
distance.matrix.pp <- vegdist(net.cut)
View(as.matrix(distance.matrix.pp))
#Mantel Test
#can do correlation, if scent similarity predicts pollinator --> mantel test
mantel(distance.matrix,distance.matrix.pp)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix, ydis = distance.matrix.pp)
Mantel statistic r: 0.02341
Significance: 0.297
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.0636 0.0844 0.1041 0.1341
Permutation: free
Number of permutations: 999
#TODO: group by family then rerun vegdist
distance.matrix.pollgroup <- vegdist(net.grp)
mantel(distance.matrix, distance.matrix.pollgroup)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix, ydis = distance.matrix.pollgroup)
Mantel statistic r: -0.06377
Significance: 0.899
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.0746 0.0954 0.1154 0.1374
Permutation: free
Number of permutations: 999
#can't predict which pollinator group will visit plant based on VOCS
#relative amounts of visits from pollinator groups
distance.matrix.pollgroup.relative <- vegdist(decostand(net.grp, method="tot"))
mantel(distance.matrix, distance.matrix.pollgroup.relative)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix, ydis = distance.matrix.pollgroup.relative)
Mantel statistic r: -0.1527
Significance: 0.985
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.107 0.140 0.170 0.217
Permutation: free
Number of permutations: 999
#similarity of the relative amounts of scent of each plant
distance.matrix.relative <- vegdist(decostand(scent.net.cut, method="tot"))
#is similarity of the relative amounts of scent of each plant correlated with the similarity in relative visits from each pollinator group
mantel(distance.matrix.relative, distance.matrix.pollgroup.relative)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix.relative, ydis = distance.matrix.pollgroup.relative)
Mantel statistic r: -0.1544
Significance: 0.985
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.110 0.143 0.172 0.205
Permutation: free
Number of permutations: 999
#No correlation between scent similarity and similarity of the proportion of visits by pollinators or pollinator group (p> 0.1)
#not grouped by pollinator class but still relative
distance.matrix.pp.relative <- vegdist(decostand(net.cut, method="tot"))
mantel(distance.matrix.relative, distance.matrix.pp.relative)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix.relative, ydis = distance.matrix.pp.relative)
Mantel statistic r: 0.05562
Significance: 0.192
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.078 0.104 0.118 0.139
Permutation: free
Number of permutations: 999
##Web plot VOCs-Plants
#interaction web plot
#plant-voc bipartite
plotweb(sqrt(scent.net), text.rot=90, method="cca", col.high=voc.class.colors[compound.class$major_class], bor.col.high=voc.class.colors[compound.class$major_class], col.low=pcols[plants.majorpoll], bor.col.low=pcols[plants.majorpoll], col.interaction=pal[5], bor.col.interaction=pal[5])
#Plant-scent-pollinator network ##NMDS
scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ]
scent.mds <- metaMDS(sqrt(scent.net.cut), autotransform = F, trace=F)
scent.op <- ordiplot(scent.mds, type="n")
points(scent.op, what="species", pch=3, col="grey50")
points(scent.op, what="sites", cex=1.5, pch=19, col=pcols[plants.majorpoll])
text(scent.op, what="sites", cex=0.8, col=pcols[plants.majorpoll])
#are scents of plants with different major pollinators different? no...
#do phylogeny, plants in same genus will cluster together need to account for phylogeny
#cap
scent.cap <- capscale(sqrt(scent.net.cut)~plants.majorpoll.cut)
anova(scent.cap)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999
Model: capscale(formula = sqrt(scent.net.cut) ~ plants.majorpoll.cut)
Df Variance F Pr(>F)
Model 3 808879 0.5422 0.792
Residual 39 19393005
plot(scent.cap, type="text")
##Indicator compounds
library(indicspecies)
plants.majorpoll.cut <- plants.majorpoll[rownames(scent.net.cut)]
mp <- multipatt(as.data.frame(scent.net.cut), plants.majorpoll.cut, control=how(nperm=999))
summary(mp)
Multilevel pattern analysis
---------------------------
Association function: IndVal.g
Significance level (alpha): 0.05
Total number of species: 355
Selected number of species: 11
Number of species associated to 1 group: 5
Number of species associated to 2 groups: 4
Number of species associated to 3 groups: 2
List of species associated to each combination:
Group butterfly #sps. 2
stat p.value
1-Hexene, 3,5-dimethyl- 0.946 0.050 *
Carbamic acid, N-phenyl-, 1,5-dimethyl-1-vinyl-4-hexenyl ester 0.942 0.038 *
Group hummingbird #sps. 3
stat p.value
Propanoic acid, 2-methyl-, 3-methylbutyl ester 1.000 0.037 *
Petasitene 0.996 0.037 *
2H-Pyran-2-one, tetrahydro-3,6-dimethyl- 0.962 0.030 *
Group butterfly+hummingbird #sps. 3
stat p.value
Bicyclo[7.2.0]undec-4-ene, 4,11,11-trimethyl-8-methylene- 0.967 0.002 **
Methoxyacetic acid, tetradecyl ester 0.941 0.015 *
1-Pentanol, 2-ethyl-4-methyl- 0.918 0.050 *
Group fly+hummingbird #sps. 1
stat p.value
2-Butenal, 3-methyl- 0.921 0.043 *
Group butterfly+fly+hummingbird #sps. 2
stat p.value
(E)-4-Oxohex-2-enal 0.941 0.012 *
1-Propanol, 2-methyl- 0.890 0.039 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Shannon diversity
###Do generalist plants have higher diversity of scent compounds?
plants.numlinks <- rowSums(decostand(net.cut, "pa"))
#polls.classtree <- compute.brlen(polls.classtree$phylo,1)
#library(picante)
#plants.pdlinks <- pd(net.lump, polls.classtree)$PD
plants.numcompounds <- rowSums(decostand(scent.net.cut, "pa"))
plants.sumcompounds <- rowSums(scent.net.cut)
plants.shannoncompounds <- diversity(scent.net.cut, index="shannon")
plants.shannonlinks <- diversity(net.grp, index="shannon")
scent.diversity <- tibble(numlinks=plants.numlinks, numcompounds=plants.numcompounds, sumcompounds=plants.sumcompounds, shannoncompounds=plants.shannoncompounds, shannonlinks=plants.shannonlinks)
ggplot(scent.diversity, aes(y=numcompounds, x=numlinks))+
geom_point()+ geom_smooth()
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
anova(lm(numlinks~numcompounds, data=scent.diversity))
Analysis of Variance Table
Response: numlinks
Df Sum Sq Mean Sq F value Pr(>F)
numcompounds 1 792.3 792.30 3.8827 0.05556 .
Residuals 41 8366.4 204.06
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#tells us if plants with more diverse scent get more pollinator species
anova(lm(plants.numlinks~shannoncompounds, data=scent.diversity))
Analysis of Variance Table
Response: plants.numlinks
Df Sum Sq Mean Sq F value Pr(>F)
shannoncompounds 1 35.4 35.438 0.1593 0.6919
Residuals 41 9123.3 222.520
#plot(plants.numcompounds~plants.pdlinks)
#plot(plants.shannoncompounds~plants.pdlinks)
shannon.compoundxlinks.plot <- ggplot(scent.diversity, aes(y=shannonlinks, x=shannoncompounds))+
geom_point()+ geom_smooth()
ggsave(shannon.compoundxlinks.plot, file='shannon.compoundxlinks.plot.jpg', width=130, height=105, units= 'mm', dpi=300)
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
compoundxlinks.plot <- ggplot(scent.diversity, aes(y=plants.numlinks, x=shannoncompounds))+
geom_point()+ geom_smooth()
ggsave(compoundxlinks.plot, file='compoundxlinks.plot.jpg', width=130, height=105, units= 'mm', dpi=300)
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
shannon.compound.plot <- ggplot(scent.diversity, aes(y=shannoncompounds, x=plants.majorpoll.cut))+
geom_point()+ geom_smooth()
ggsave(shannon.compound.plot, file='shannon.compound.plot.jpg', width=110, height=105, units= 'mm', dpi=300)
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
numcompounds.plot <-ggplot(scent.diversity, aes(y=numcompounds, x=plants.majorpoll.cut))+
geom_point()+ geom_smooth()
ggsave(numcompounds.plot, file='numcompounds.plot.jpg', width=110, height=105, units= 'mm', dpi=300)
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
sumcompounds.plot<- ggplot(scent.diversity, aes(y=sumcompounds, x=plants.majorpoll.cut))+
geom_point()+ geom_smooth()
ggsave(sumcompounds.plot, file='sumcompounds.plot.jpg', width=110, height=105, units= 'mm', dpi=300)
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
#VOC-pollinator network
#interaction web plot
#voc-pollinator bipartite
#make dataframe with vocs and pollinators
##change scent.net to long format
plantcompoundemission <- scent.net.cut %>% as.data.frame() %>% rownames_to_column("plant") %>% #average peak area
pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)
#write_csv(plantcompoundemission, "plantcompoundemissions.csv")
#total number of visits to plant by each pollinator species
plantpollinatorvisits <- net.cut %>% as.data.frame() %>% rownames_to_column("plant") %>% pivot_longer(-plant, names_to = "pollinator", values_to = "visits")
#total number of visits to each compound. Pollinators rows, compounds cols, values = number of visits
scentpollnet <- left_join(plantcompoundemission, plantpollinatorvisits,relationship = "many-to-many") %>% group_by(pollinator, compound) %>% summarise(visits=sum(visits)) %>% pivot_wider(names_from = "compound", values_from = "visits") %>% column_to_rownames("pollinator") %>% as.matrix()
Joining with `by = join_by(plant)`
`summarise()` has grouped output by 'pollinator'. You can override using the
`.groups` argument.
#relative amount of each compound
plantcompoundrelative <- scent.net.cut %>% as.data.frame() %>% decostand(method="tot") %>% rownames_to_column("plant") %>% #average peak area
pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)
#multiply number of visits by relative peak area
# weights visit to compounds. Proportions visits based on how many VOCs there are.
scentpollnet.area <- left_join(plantcompoundrelative, plantpollinatorvisits,relationship = "many-to-many") %>% group_by(pollinator, compound) %>% mutate(visits.area=visits*emissions) %>% summarise(visits.area=sum(visits.area)) %>% pivot_wider(names_from = "compound", values_from = "visits.area") %>% column_to_rownames("pollinator") %>% as.matrix()
Joining with `by = join_by(plant)`
`summarise()` has grouped output by 'pollinator'. You can override using the
`.groups` argument.
scentpollnet.class <- tibble(pollinator=rownames(scentpollnet)) %>% left_join(polls.class.dataframe) %>% pull(pollinator_group)
Joining with `by = join_by(pollinator)`
#plot bipartite network
#plotweb(sqrt(scentpollnet), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], #col.low=pcols[scentpollnet.class], bor.col.low=pcols[scentpollnet.class], col.interaction=pal[5], #bor.col.interaction=pal[5])
#plotweb for scentpollnet.area
plotweb(sqrt(scentpollnet.area), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], col.low=pcols[scentpollnet.class], bor.col.low=pcols[scentpollnet.class], col.interaction=pal[5], bor.col.interaction=pal[5])
##Heatmap scent-net
#value = number of
heatmaply(scentpollnet, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
#row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
#row_side_palette=function(n) pcols,
row_side_colors=data.frame(PollinatorGroup=scentpollnet.class),
row_side_palette=function(n) pcols.all)
voc.pollgrp.cap <- capscale(sqrt(scentpollnet.area)~scentpollnet.class)
anova(voc.pollgrp.cap)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999
Model: capscale(formula = sqrt(scentpollnet.area) ~ scentpollnet.class)
Df Variance F Pr(>F)
Model 4 15.904 1.67 0.105
Residual 88 209.505
plot(voc.pollgrp.cap, type="text")
#chemical class x pollinator group bipartite
#make vector with compound class and all scent compounds
vol.class <-deframe(select(compound.class, compound_long, major_class) )
#turns scentpollnet matrix to data frame then pivots longer to have 2 cols called chemical and visits
chemical.taxa <- left_join(plantcompoundemission, plantpollinatorvisits,relationship = "many-to-many") %>% left_join(polls.class.dataframe) %>% left_join(compound.class, by=c("compound"= "compound_long")) %>% group_by(pollinator_group, major_class) %>% summarise(visits=sum(visits)) %>% pivot_wider(names_from = major_class, values_from = visits) %>% column_to_rownames("pollinator_group") %>% as.matrix()
Joining with `by = join_by(plant)`
Joining with `by = join_by(pollinator)`
`summarise()` has grouped output by 'pollinator_group'. You can override using
the `.groups` argument.
#chemical class & pollinator group network
# how do I incorporate chemical class instead of compounds?
#make data frame with chemical class and pollinator family - use scent.net & chemical.class
#and use scentpollnet but scentpollnet compounds in different order than scent.net
plotweb(sqrt(chemical.taxa), text.rot=90, method="normal", col.high=pal[3], bor.col.high=pal[3], col.low=pcols.all, bor.col.low=pcols.all, col.interaction=pal[5], bor.col.interaction=pal[5])
library(bipartite)
library(vcd)
Loading required package: grid
#side quest
mosaic(chemical.taxa)
mosaic(net)
chisq.test(chemical.taxa) #p=1 so this means that no specialist compound classes. and vise versa. nobody is specialized
Pearson's Chi-squared test
data: chemical.taxa
X-squared = 8488.8, df = 24, p-value < 2.2e-16
#relative peak area x visits
#value = number of
heatmaply(scentpollnet.area, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
#row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
#row_side_palette=function(n) pcols,
row_side_colors=data.frame(PollinatorGroup=scentpollnet.class),
row_side_palette=function(n) pcols.all)
#Rarifaction
#species accumulation curves
#figure out what percent Pauls data catpured the whole pollinator community
#pull out 1 random 1 hr sampling period and say how many pollinators observered, then do this again and again randomly - variation in samples = standard error y axis= number of species, x= number of pollinator observations
#class= bee, fly, butterfly, hummingbird
#xaxis is number of plant species observed (if you did the whole summer but only looked at 1 species and then next summer looked at 2 species, etc.)
#
#simplify net.long to only have numbers in it so that accumulation curve will work
#rows are observation "units" and columns are pollinators
net.wide <- net.long %>% select(pollinator, interactions) %>% mutate(index=row_number()) %>%
pivot_wider(names_from=pollinator, values_from = interactions) %>% mutate(across(everything(),~replace_na(.x, 0))) %>% select(-index) %>% as.matrix()
for(class in unique(polls.class)) { #could replace class with family
class_subset <- net[,polls.class==class] #making accumulation curve for each pollinator class
class_sa <- vegan::specaccum(class_subset, permutations = 999)
plot.new()
plot(class_sa, col = as.integer(factor(class))+1, add = class != "bee", lwd=2,
xlab="Pollinator Observations", ylab="Pollinator Species", main="Pollinator Species accumulation curve", sub=class)
}
legend("bottomright", unique(polls.class),
col = 2:5, pch = 19, title = "Class")
sa<- vegan::specaccum(net.wide,permutations=100, method="exact")
plot(sa)
#Pauls total pollinator observations 8hrs/wk * 42wks = 336 hours total
visitationrate <- nrow(net.long)/336 #15.33 unique p-p observations per hr (observation can include multiple visits)
#observation=within a 15min period how many p-p interactions did we see, total number of links observed
2000/visitationrate # 130 hrs per treatment ( 2 treatments = 260 hrs per year)
[1] 130.4854
#then based on my observation hours, what percent of the pollinator community could be captured
#normalized degree for p-p
library(rstatix)
Attaching package: 'rstatix'
The following object is masked from 'package:stats':
filter
ND_poll<-ND(net, normalised = T)$higher
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
ND.poll.class <- tibble(ND=ND_poll, polls.class=polls.class)
kruskal_test(ND.poll.class, ND~polls.class)
# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <int> <dbl> <chr>
1 ND 93 4.01 4 0.404 Kruskal-Wallis
ggplot(ND.poll.class, aes(x= polls.class, y=ND)) +geom_boxplot()
#normalized degree for VOC-pollinator
library(rstatix)
#ND for pollinators in scent network
ND_poll_scent<-ND(scentpollnet, normalised = T)$lower
#ND of VOC chemical class
ND_scent_poll<-ND(scentpollnet, normalised = T)$higher
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
ND.poll <- tibble(ND=ND_poll_scent, polls.class=polls.class)
kruskal_test(ND.poll, ND~polls.class)
# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <int> <dbl> <chr>
1 ND 93 2.69 4 0.611 Kruskal-Wallis
ggplot(ND.poll, aes(x= polls.class, y=ND)) +geom_boxplot()
#creates dataframe with normalized degree for each VOC and chemical class (same as data.frame)
ND.scent <- tibble(ND=ND_scent_poll, vols.class=compound.class$major_class)
#as.data.frame(list(ND = ND_scent_poll, vols.class = compound.class$major_class))
kruskal_test(ND.scent, ND~vols.class)
# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <int> <dbl> <chr>
1 ND 354 7.83 6 0.251 Kruskal-Wallis
ggplot(ND.scent, aes(x= vols.class, y=ND)) +geom_boxplot()
#normalized visits - home cooked
#created a function similar to normalized degree(ND) but ND looks at presence/absence in the network where normalized visits (NV) looks at the normalized weighted visits/interactions. This accounts for variation in the amount each plant or compound is visited and give those with higher interaction frequency greater weight or contribution to the network
normalized.visits <- function (web, normalised = TRUE)
{
websum <- sum(web) #takes sum of all interactions in the web
dlower <- rowSums(web)
dhigher <- colSums(web)
Nlow <- Nhigh <- 2
if (normalised) {
Nlow <- websum
Nhigh <- websum
}
low <- dlower/Nlow #weighting the interactions by dividing rowsums by the sum of all links in the network
names(low) <- rownames(web) #dividing colsums by the sum of all links in the network
high <- dhigher/Nhigh
names(high) <- colnames(web)
list(lower = low, higher = high)
}
#weighted normalized visits for p-p (pollinators only)- home cooked
#Which pollinators have the most visits. of all the visits recorded, what percentage of them were by that pollinator (for a given plant what percentage of visits does it get out of the whole meadow, for a pollinator what percentage of visits does it give out of all visits)
#normalized visits for pollinators in scent net
NV_poll<-normalized.visits(net, normalised = T)$higher
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.poll.class <- tibble(NV=NV_poll, polls.class=polls.class)
kruskal_test(NV.poll.class, NV~polls.class)
# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <int> <dbl> <chr>
1 NV 93 11.6 4 0.0208 Kruskal-Wallis
ggplot(NV.poll.class, aes(x= polls.class, y=NV)) +geom_boxplot()
#only a few species making most of the visits
#1 bee responsible for over 40% of visits
#normalized visits for scent network
#normalized visits for pollinators in scent net
#which pollinator class visits the most amount of compounds weighted by the number of visits
NV_poll_scentnet<-normalized.visits(scentpollnet, normalised = T)$lower
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.poll.class.scentnet <- tibble(NV=NV_poll_scentnet, polls.class=polls.class)
kruskal_test(NV.poll.class.scentnet, NV~polls.class)
# A tibble: 1 × 6
.y. n statistic df p method
* <chr> <int> <dbl> <int> <dbl> <chr>
1 NV 93 2.91 4 0.573 Kruskal-Wallis
ggplot(NV.poll.class.scentnet, aes(x= polls.class, y=NV)) +geom_boxplot()
#super generalist bees
#normalized visits for VOCs in scent net
#which compound class are most visited
NV_voc<-normalized.visits(scentpollnet, normalised = T)$higher
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.vols.class <- tibble(NV=NV_voc, vols.class=compound.class$major_class)
#TODO: Fix
k_t.vols.class <- kruskal_test(NV.vols.class, NV~vols.class)
library(coin)
Loading required package: survival
Attaching package: 'coin'
The following objects are masked from 'package:rstatix':
chisq_test, friedman_test, kruskal_test, sign_test, wilcox_test
The following object is masked from 'package:sna':
rperm
library(PMCMRplus) #for posthoc
library(PMCMR)
PMCMR is superseded by PMCMRplus and will be no longer maintained. You may wish to install PMCMRplus instead.
#posthoc_vols.class <- posthoc.kruskal.nemenyi.test(NV.vols.class$vols.class, NV.vols.class$NV, dist = "tukey")
ggplot(NV.vols.class, aes(x= vols.class, y=NV)) +geom_boxplot()
#tells you what vocs get the most visits
#no difference in visitation rate to compound classes
#modularity for net
library(vcd)
modules.net <- computeModules(sqrt(net))
plotModuleWeb(modules.net, labsize=.5)
module.net.list <- listModuleInformation(modules.net)
flatten(module.net.list[[2]])
[[1]]
[1] "androsace_septentrionalis" "claytonia_lanceolata"
[3] "draba_aurea" "erythronium_grandiflorum"
[5] "hydrophyllum_capitatum" "lomatium_dissectum"
[7] "pseudocymopterus_montanus" "ranunculus_inamoenus"
[9] "taraxacum_officinale" "viola_praemorsa"
[11] "amelanchier_alnifolia" "fragaria_virginiana"
[13] "hydrophyllum_fendleri" "linum_lewisii"
[[2]]
[1] "andrena_cyanophila" "lasioglossum_spp" "sphecodes_sp_01"
[4] "syrphidae_spp" "scathophagidae_sp_01" "halictus_confusus"
[7] "halictus_virgatellus" "halictidae_green_spp_1" "eupeodes_volucris"
[10] "andrena_sp_02" "muscidae_sp_02" "syrphidae_02"
[13] "rhamphomyia_sp_01" "cartosyrphus_tarda" "callophrys_affinis"
[16] "osmia_lignaria" "cartosyrphus_sp_02" "tephritidae_01"
[19] "polygonia_zephyrus"
[[3]]
[1] "delphinium_nuttallianum" "ipomopsis_aggregata"
[3] "gentiana_parryi" "vicia_americana"
[[4]]
[1] "bombus_appositus" "selasphorus_platycercus"
[3] "anthidium_emarginatum" "bombus_nevadensis"
[5] "papilio_gothica"
[[5]]
[1] "agoseris_glauca" "arenaria_congesta" "eriogonum_subalpinum"
[4] "eriogonum_umbellatum" "galium_boreale" "mahonia_repens"
[7] "potentilla_hippiana" "rosa_woodsii" "sedum_lanceolatum"
[10] "achillea_millefolium" "potentilla_gracilis"
[[6]]
[1] "muscidae_spp_01" "colletes_kincaidii"
[3] "halictus_rubicundus" "halictidae_spp_01"
[5] "panurginus_ineptus" "thricops_septentrionalis"
[7] "osmia_bucephala" "chrysis_coerulans"
[9] "villa_eumenes" "sphaerophoria_robusta"
[11] "glaucopsyche_lygdamus" "pieris_mcdunnoughi"
[13] "andrena_vicinoides" "colletes_nigrifrons"
[15] "melanstoma_caerulescens"
[[7]]
[1] "agoseris_aurantiaca" "campanula_rotundifolia"
[3] "erigeron_speciosus" "helianthella_quinquenervis"
[5] "heliomeris_multiflora" "heterotheca_villosa"
[7] "senecio_serra" "solidago_multiradiata"
[9] "castilleja_sulphurea" "phacelia_heterophylla"
[[8]]
[1] "nymphalidae_spp" "bombus_bifarius"
[3] "megachile_relativa" "tachinidae_sp_01"
[5] "bombus_sylvicola" "lycaenidae_spp"
[7] "anthophora_terminalis" "bombus_flavifrons"
[9] "andrena_costillensis" "bombus_californicus"
[11] "bombus_occidentalis" "megachile_pugnata"
[13] "arctophila_flagrans" "coelioxus_funeraria"
[15] "osmia_grindeliae" "pieridae_spp"
[17] "speyeria_mormonia" "chrysotoxum_ventricosum"
[19] "megachile_melanophea" "psithyrus_insularis"
[21] "lycaena_rubidus_sirius" "colias_alexandra"
[23] "bombus_frigidus" "mesebrina_latreillei"
[25] "syrphidae_spp_2" "syrphidae_spp_4"
[27] "sphex_spp_1" "bombus_balteatus"
[29] "syrphidae_spp_1" "phormia_spp_1"
[31] "megachile_frigida" "lycaena_cupreus"
[33] "bombus_mixtus" "osmia_tristela"
[35] "melissodes_confusa"
[[9]]
[1] "boechera_stricta" "erigeron_flagellaris" "lathyrus_lanszwertii"
[4] "mertensia_ciliata" "mertensia_fusiformis" "senecio_integerrimus"
[[10]]
[1] "ochlodes_sylvanoides" "systoechus_sp" "osmia_proxima"
[4] "pieris_rapa" "andrena_transnigra" "hoplitus_fulgida"
[7] "osmia_iridis" "osmia_subaustralis" "osmia_montana"
[10] "osmia_coloradensis" "melanostoma_kelloggi" "syrphidae_spp_3"
[13] "eupeodes_lapponicus" "euphydras_anicia" "eristalis_latifrons"
[16] "osmia_sp_sgo" "cryptopogon_sp_01" "hoplitus_robusta"
[19] "coenonympha_ochracea"
pollinator.modules <- tibble(module=1:5) %>% mutate(pollinator=map(module, ~module.net.list[[2]][[.x]][[2]])) %>% unnest(pollinator) %>% left_join(polls.class.dataframe) %>% count(module, pollinator_group) %>% pivot_wider(names_from = pollinator_group, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()
Joining with `by = join_by(pollinator)`
chisqr.polls <- chisq.test(pollinator.modules)
Warning in chisq.test(pollinator.modules): Chi-squared approximation may be
incorrect
residules.chisqr.polls <- residuals(chisqr.polls)
mosaic(t(pollinator.modules), shade = T, gp=shading_hcl)
#test to see p-p network modules
#mosaic(net, shade = T, gp=shading_hcl)
#modularity for scent net
modules.scentnet <- computeModules(sqrt(t(scentpollnet)))
#plotModuleWeb(modules.scentnet)
module.scentnet.list <- listModuleInformation(modules.scentnet)
pollinator.modules <- tibble(module=1:4) %>% mutate(pollinator=map(module, ~module.scentnet.list[[2]][[.x]][[2]])) %>% unnest(pollinator) %>% left_join(polls.class.dataframe) %>% count(module, pollinator_group) %>% pivot_wider(names_from = pollinator_group, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()
Joining with `by = join_by(pollinator)`
chisqr.polls <- chisq.test(pollinator.modules)
Warning in chisq.test(pollinator.modules): Chi-squared approximation may be
incorrect
residules.chisqr.polls <- residuals(chisqr.polls)
mosaic(t(pollinator.modules), shade = T, gp=shading_hcl, labeling = vcd::labeling_border(rot_labels = c(0, 0)))
#chemical class modularity
chemical.modules <- tibble(module=1:4) %>% mutate(compound_long=map(module, ~module.scentnet.list[[2]][[.x]][[1]])) %>% unnest(compound_long) %>% left_join(compound.class) %>% count(module, major_class) %>% pivot_wider(names_from = major_class, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()
Joining with `by = join_by(compound_long)`
mosaic(t(chemical.modules), shade = T, gp=shading_hcl, labeling = vcd::labeling_border(rot_labels = c(0, 0)))
#stacked bar chart with flowers visited by each pollinator
library(ggplot2)
library(tidyverse)
library(RColorBrewer)
mypalette= c(brewer.pal(8, "Set2"), brewer.pal(12, "Set3"))
#looks at top 20 plants and top 75 pollinators (interactions sqrt to adjust scale on figure)
#pollinator perspective
ggplot(net.long %>% mutate(plant=fct_lump_n(plant, 18), pollinator=fct_lump_n(pollinator, 75)), aes(x=sqrt(interactions), y=pollinator, fill=plant)) + geom_col() + scale_fill_manual(values=mypalette)
#plant perspective
ggplot(net.long %>% mutate(plant=fct_lump_n(plant, 45), pollinator=fct_lump_n(pollinator, 18)), aes(x=sqrt(interactions), y=plant, fill=pollinator)) + geom_col() + scale_fill_manual(values=mypalette)
#phenology
phenology <- read_sheet("15Jnx8OWMzpakC0YvJWmZMmaRSnA0X-7HBjymDsvUfjY")
! Using an auto-discovered, cached token.
To suppress this message, modify your code or options to clearly consent to
the use of a cached token.
See gargle's "Non-interactive auth" vignette for more details:
<https://gargle.r-lib.org/articles/non-interactive-auth.html>
ℹ The googlesheets4 package is using a cached token for 'jbohey@uci.edu'.
✔ Reading from "caradonna_rmbl_flowering_phenology_data_EDI".
✔ Range 'caradonna_rmbl_flowering_phenology_data_EDI'.
#stack plots and heatmaps
#Flowers#
ggplot(phenology, aes(x=week_num, y=flower_count, fill= plant )) + geom_col() + scale_fill_manual(values=sample(rainbow(46))) +facet_wrap(vars(year))
#heatmap
ggplot(phenology, aes(x=week_num, y=plant, fill= sqrt(flower_count ))) + geom_tile() + scale_fill_viridis_c() + facet_wrap(vars(year))
#pollinators#
ggplot(net.long, aes(x=week_num, y=interactions, fill= pollinator )) + geom_col() + scale_fill_manual(values=sample(rainbow(93))) +facet_wrap(vars(year))
#heatmap
ggplot(net.long, aes(x=week_num, y=pollinator, fill= sqrt(interactions ))) + geom_tile() + scale_fill_viridis_c() + facet_wrap(vars(year))
#scent#
#number of flowers per week and multiple by mean scent of that species
#how many flowers are open each week
flowers.per.week <- phenology %>% group_by(year, week_num, plant) %>%
summarise(flower_count=sum(flower_count))
`summarise()` has grouped output by 'year', 'week_num'. You can override using
the `.groups` argument.
#ideally want for each week, 1 average weighted mean of scent compound
#but we only have a few scent samples so just taking the mean of scent compounds for each plant
#flowers per week percentage
flower.percentage <- flowers.per.week %>% group_by(year, week_num) %>%
mutate(flower_count=flower_count/sum(flower_count)) %>%
pivot_wider(names_from = c(year,week_num), values_from = flower_count, values_fill = 0) %>%
filter(plant%in%rownames(scent.net)) %>%
arrange(match(plant, rownames(scent.net)))%>% #matches scent.net plant order
column_to_rownames("plant") %>% as.matrix()
#change names in pauls data
scentxphenology <- t(flower.percentage) %*% scent.net[rownames(flower.percentage),]
#values in matrix = weighted average of flower scents and mean is weighted by how many flowers there are that week
#replaces Nans with 0
scentxphenology[is.nan(scentxphenology)] <- 0
#ggplot likes long format
scent.phenology <- scentxphenology %>% as.data.frame() %>% rownames_to_column("yearweek") %>% separate(yearweek, into = c("year", "week_num")) %>% mutate(year=as.integer(year), week_num= as.integer(week_num)) %>% pivot_longer(all_of(colnames(scent.net)))
ggplot(scent.phenology, aes(x=week_num, y=value, fill= name )) + geom_col() + scale_fill_manual(values=sample(rainbow(355)), guide="none") +facet_wrap(vars(year))
#heatmap
ggplot(scent.phenology, aes(x=week_num, y=name, fill= sqrt(value ))) + geom_tile() + scale_fill_viridis_c() + facet_wrap(vars(year))
#mega table ##flower and poll traits ###VOCS
#number of compounds each plant species has
#total emissions for each species
#diversity of scents (shannon diversity)
library(tidyverse)
scent.net.long <- scent.net %>% as.data.frame() %>% rownames_to_column("plant") %>% #average peak area
pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)
compound.class.percentage <- scent.net.long %>% left_join(compound.class, by=c(compound="compound_long")) %>% group_by(plant, major_class) %>% summarise(emissions=sum(emissions)) %>% #gives total emissions for each compound class
mutate(emissions=emissions/sum(emissions)) %>% #relative amount of each compound class in species scent
pivot_wider(names_from = major_class, values_from = emissions, names_prefix = "percent_", values_fill = 0) #makes column for each compound class and gives percentage of that compound class out of the total emissions for that species
`summarise()` has grouped output by 'plant'. You can override using the
`.groups` argument.
#rare vs common compounds (mean rarity)
##for each compound, how many times does it occur in the dataset?
compound.rarity <- scent.net.long %>% count(compound) %>% #number of species the compound occurs in
mutate(rarity=n/nrow(scent.net)) #the percentage of species the compound occurs in
weighted.scent.rarity <- scent.net.long %>% left_join(compound.rarity) %>% group_by(plant) %>%
mutate(emissions=emissions/sum(emissions)) %>% #makes emissions -> relative emissions
mutate(weighted_scent_rarity = emissions*rarity) %>% #multiplies emissions x rarity to get weighted average
summarise(weighted_scent_rarity= sum(weighted_scent_rarity)) #then add the different parts of the average
Joining with `by = join_by(compound)`
#rarity score of 1 means this compound is found in all the plant species, 0= no other species has these compounds
#mertensia has 0.70 rarity - the weighted average of all the compounds rarity; 70% of compounds (or total emissions) are shared, 30% of compounds or total emissions other species don't have
scent.summary <- tibble(plant=rownames(scent.net), total_emissions= rowSums(scent.net), number_compounds=rowSums(scent.net>0),scent_shannon_diversity= diversity(scent.net, index="shannon")) %>%
left_join(weighted.scent.rarity) %>%
left_join(compound.class.percentage)
Joining with `by = join_by(plant)`
Joining with `by = join_by(plant)`
##floral traits summary
#pollinators that visit plants
#Percentage of visits by each pollinator family to each plant species
pollinator.family.percentage <- net.long %>% group_by(plant, pollinator_family) %>% summarise(interactions=sum(interactions)) %>% #gives sum of interactions for each pollinator family
mutate(interactions=interactions/sum(interactions)) %>% #relative amount of each pollinator family visiting each flower species
pivot_wider(names_from = pollinator_family, values_from = interactions, names_prefix = "percent_", values_fill = 0)
`summarise()` has grouped output by 'plant'. You can override using the
`.groups` argument.
#Percentage of visits by each pollinator group to each plant species
pollinator.group.percentage <- net.long %>% group_by(plant, pollinator_group) %>% summarise(interactions=sum(interactions)) %>% #gives sum of interactions for each pollinator group
mutate(interactions=interactions/sum(interactions)) %>% #relative amount of each pollinator group visiting each flower species
pivot_wider(names_from = pollinator_group, values_from = interactions, names_prefix = "percent_", values_fill = 0)
`summarise()` has grouped output by 'plant'. You can override using the
`.groups` argument.
#rarity - how generalized the pollinator is
pollinator.rarity <- net.long %>% count(plant, pollinator) %>% #adds how many rows of interactions between plant & poll species
count(pollinator, name="number_plants_visited") %>% #number of plants the pollinator goes to
mutate(rarity=number_plants_visited/nrow(net)) #the percentage of plant species the pollinator visits in the meadow
#gives for each plant species, which pollinators visited and their relative contributions to total number of visits
net.long.relative.interactions <- net.long %>% group_by(plant, pollinator) %>%
summarise(interactions=sum(interactions)) %>%
left_join(pollinator.rarity) %>%
group_by(plant) %>%
mutate(interactions=interactions/sum(interactions))
`summarise()` has grouped output by 'plant'. You can override using the
`.groups` argument.
Joining with `by = join_by(pollinator)`
#percent of interactions that come from each pollinator species
weighted.pollinator.rarity <- net.long.relative.interactions %>% mutate(weighted_poll_rarity = interactions*rarity) %>% #multiplies interactions x rarity to get weighted average
summarise(weighted_poll_rarity = sum(weighted_poll_rarity))
flower.visitor.summary <- tibble(plant=rownames(net), total_visits= rowSums(net), number_pollspecies_visited=rowSums(net>0),pollinator_shannon_diversity= diversity(net, index="shannon")) %>%
left_join(weighted.pollinator.rarity) %>%
left_join(pollinator.family.percentage) %>% left_join(pollinator.group.percentage)
Joining with `by = join_by(plant)`
Joining with `by = join_by(plant)`
Joining with `by = join_by(plant)`
megatable <- scent.summary %>% full_join(flower.visitor.summary) #full join bc pauls data and scent data plants don't match
Joining with `by = join_by(plant)`
#correlations
library(pheatmap)
library(smatr)
correlation.table <- megatable %>% column_to_rownames("plant") %>% cor(use="na.or.complete")
pheatmap(correlation.table, scale="none", clustering_method = "mcquitty", color=colorRampPalette(c("blue","white", "red"))(200),
breaks = seq(-1,1, by=0.01)) #which ones it labels on the scale
#scatterplot with just weighted scent and weighted pollinator visits
weighted.scent.poll <- ggplot(megatable, aes(x=weighted_poll_rarity, y=weighted_scent_rarity, color=number_pollspecies_visited)) + geom_point()+geom_smooth() + geom_text(aes(label=plant))
print(weighted.scent.poll)
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Warning: Removed 5 rows containing non-finite outside the scale range
(`stat_smooth()`).
Warning: The following aesthetics were dropped during statistical transformation:
colour.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_point()`).
Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_text()`).
#the most common volatiles get the most visits from the most generalized pollinators
sma.weighted.scent.poll <- sma(weighted_poll_rarity~weighted_scent_rarity, data=megatable )
summary(sma.weighted.scent.poll)
Call: sma(formula = weighted_poll_rarity ~ weighted_scent_rarity, data = megatable)
Fit using Standardized Major Axis
------------------------------------------------------------
Coefficients:
elevation slope
estimate -0.3602556 1.352232
lower limit -0.5885207 1.001742
upper limit -0.1319906 1.825352
H0 : variables uncorrelated
R-squared : 0.06769348
P-value : 0.091989
#percent_hesperiidae x nitrogenous vocs
ggplot(megatable, aes(x=percent_nitrogenous, y=percent_lycaenidae, color=number_pollspecies_visited)) + geom_point()+geom_smooth() + geom_text(aes(label=plant))
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Warning: Removed 5 rows containing non-finite outside the scale range
(`stat_smooth()`).
Warning: The following aesthetics were dropped during statistical transformation:
colour.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_point()`).
Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_text()`).
##pollinator trait summary